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Abstract. As an introductory lecture to the workshop an overview is given 
over continuum models for homoepitaxial surface growth using partial dif- 
ferential equations (PDEs). Their heuristic derivation makes use of inherent 
symmetries in the physical process (mass conservation, crystal symmetry, . . . ) 
which determines their structure. Two examples of applications are given, one 
for large scale properties, one including crystal lattice discreteness. These are: 
(i) a simplified model for mound coarsening and (ii) for the transition from 
layer-by-layer to rough growth. Virtues and shortcomings of this approach is 
discussed in a concluding section. 
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1. Introduction 

Crystal growth by Molecular Beam Epitaxy (MBE) involves processes on quite 
different physical scales, from atomic to micrometer scales. To capture the large 
scale features it may be appropriate to represent the surface as a continuous height 
field /i(x, t) above a set of so called base points x where h : (x, t) 6 M 2 x R — > R. 
The surface dynamics is then modelled by a parabolic PDE for ft,(x, i), which 
often is a starting point interesting mathematical problems, and in some cases for 
rigorous results 8,9, which were an important part of this workshop. 

Their mathematical tractability somehow contrasts with the lack of rigour 
in their derivations, as we shall see below. It would be nice to "prove" a certain 
PDE to be the appropriate large scale model for some underlying more elementary 
model, typically a stochastic model of the solid-on-solid type. This would connect 
the result for large scales with the microscopic model which is much closer to 
physical reality. 

Nevertheless, often this chain can be closed in a heuristic sense, and this 
lecture presents some examples. It is organised as follows: First, the heuristic pro- 
cedure of deriving a continuum model is motivated and its structure is explained. 
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Then, two examples of applications are presented, covering in some sense extreme 
cases of applicability of continuum models; (i) coarsening of mound patterns in 
homoepitaxy and (ii) damping of width oscillations in the transition from layer- 
by-layer to rough growth. Last, the general applicability of continuum models to 
MBE growth is discussed, in particular those features of discreteness of the un- 
derlying process which can only be covered incompletely by a continuum model. 
Open questions in the derivation which have not been addressed so far, but appear 
to be solvable, are mentioned here. 

The main purpose of this introductory lecture is to recollect the physical 
motivation of the class of models which deal with the largest spatial and tem- 
poral scales involved in MBE growth, and to give some intuition to other more 
mathematical approaches. 

2. Derivation and structure 
2.1. General background 

Continuum models for surface growth are comparable to the hydrodynamic limit of 
microscopic models of fluid dynamics. Behind the continuous macroscopic height 
field /i(x, t) one has to imagine an ensemble of microscopic configurations which 
are locally in a steady state. Figure ^ illustrates this relationship. The local state 
is subject to changes on time scales, ideally longer than the relaxation times into 
those states, which leads to the dynamics of the height field. 

The structure of the PDE governing the dynamics of h(x, t) reflects differ- 
ent elementary processes and their symmetries. In growth from molecular beam 
epitaxy there is 

• a deposition flux F^y (x, t), i.e., the deposited volume per time and surface 
area, 

• possibly an evaporation loss E^ h y(x, t), again volume per time and area, 

• a surface current J{^-j(x, t) on terraces and along steps, transported volume 
per time and length, i.e., cross section on the surface. 

Multiplied with a volume density p these quantities become mass fluxes and cur- 
rents. All three depend on the actual surface configuration h, e.g., on local slopes, 
their orientations and curvatures. Elastic interactions and terrace adatom detach- 
ment with subsequent redeposition couple the dynamics at a surface point non- 
locally to surface configurations at other places. There may also be an explicit 
space dependence beyond the influence of the height configuration, e.g., for a non- 
uniform deposition with pulsed flux. 

Deposition, evaporation, and surface current are combined into a continuity 
equation 

— / i (x,t)=f {fi) (x,()-V.J w (x,t)-%(x,t). (2.1) 

Note that small local slopes are assumed, such that any geometric prefactors to 
the divergence of the surface current typically are neglected. 
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Figure 1 . Sketch of the hight field representation for surface con- 
figurations. Upper: one surface configuration out of an ensemble 
in local steady states. Lower: a smooth height field h(x, t) as a 
representation. 



A full derivation from local steady states in the sense postulated above is of 
course impossible in almost any case, apart from some one-dimensional examples 
P). One therefore has to rely on heuristically justified (reasonable) assumptions 
for the dynamic quantities involved. The connection to the more detailed models 
on smaller scales is given by length scales which influence the surface morphology 
and become visible. Explained more precisely in the other introductory lectures of 
this conference *^ ey are 

• a, the lattice spacing of the crystal. 

• £d, the diffusion length, i.e., the "typical" distance of nucleation sites on a 
terrace, with a non-trivial dependence on the deposition rate F and adatom 
diffusivity D. 

• £~es, the Ehrlich-Schwoebel length, i.e., a measure for the reflectivity of the 
barrier at a downward step. Roughly speaking, an adatom located right at a 
downward step is equally likely to hop downward as it is to reach the next 
upward step, if it lies at a distance £es- 

• ^step, the average distance between two truly ascending at a vicinal surface 
of slope a/£ s tep without counting islands between them. 
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2.2. Surface current 

Ehrlich Schwoebel effect on terraces 

As was realised quite early the most important ingredient of Eq. 1|2.1[) is a non- 
equilibrium current on surface terraces induced by Ehrlich-Schwoebel (ES) barriers 

EH EH EH 

For small average slopes m = V/i with |m| <C o/£d or £b -C 4tep ; the 
ES barrier induces a net uphill current Jest ~ Fm. Inserted into the surface 
current in Eq. (|2.1(l this results in a destabilising term for the height equation 
d t h = -cV -Vh+ .... 

Of course one would like to know the functional dependence of the surface 
current Jest beyond this approximation for small slopes. Following an argument 
of Siegert for 1 + 1 dimensions \Fij\ there must also be surface orientations with 
stabilising surface current: Think of Jest as a vector field on the tangent bundle 
of the unit sphere of possible surface orientations. By the ES instability all high 
symmetry orientations have Jest pointing outward around a zero. By continuity 
of Jest as a function of the orientation there must be zeroes with the current field 
pointing inward to the zero, i.e., orientations stabilised by the surface current. 

These qualitatively fundamental properties, the existence of unstable and 
stable surface orientations, can be cast into a heuristically postulated function 
Jest(hi). Figure |2 shows two such examples, one isotropic within the crystal 




Figure 2. Examples of heuristic choices for nonequilibrium ter- 
race Ehrlich-Schwoebel surface current Jest(hi)- Left: in-plane 
isotropy, Jest(hi) = (1 — l m | 2 ) m : continuous set of stable 
slopes {m | |m| = 1}. Right: anisotropy, four-fold symmetry, 
JEST,i(m) = (1 — mf) vrii with i = x,y, discrete set of stable 
slopes at (±1, ±1). 

terrace — a somewhat artificial but nevertheless theoretically interesting case, one 
isotropic, with four-fold (7r/2) rotational symmetry. 
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Ehrlich Schwoebel effect at steps 

The Ehrlich-Schwoebel effect also has an influence on adatom diffusion along 
steps, allowing for attachment of atoms to kinks preferentially from one side, as 
is sketched in Fig. |3J ft is one important mechanism causing a step meandering 



Figure 3. The ES effect at step kinks favours adatom incorpo- 
ration from one side and induces a net current along the step. 

instability on growing vicinal surfaces jf 31 [7] . Politi and Krug derived the contri- 
bution Jess of non-equilibrium current along steps to the total surface current for 
a surface with given average slope which is assumed to be constant on a sufficiently 
large region such that steps have a well defined average density and orientation 
with respect to the underlying crystal lattice |f 5| . 

Curvature relaxation 

Clearly on small scales the destabilising ES effects encountered in this section so 
far are balanced. Too thin protrusions are forbidden by their cost in surface free 
energy. For thermal relaxation towards equilibrium the corresponding continuum 
dynamics are straightforward to derive. The surface free energy 

T= J j(Vh) y/1 + (Vh) 2 d d x (2.2) 

weights each surface element by its orientation dependent contribution -f(Vh). The 
gradient in chemical potential /i = SJ- / 5h drives the relaxational surface current 
which is also proportional to the (tensorial and orientation dependent) adatom 
mobility <j(Vh). A linear expansion around a surface of constant slope h = m • x 
small enough to again neglect geometric effects |2H| yields 

Jcurv = o- V V-[ 7 1 + (V m V m 7)]-V/i. (2.3) 

The expression in square brackets is generally called the surface stiffness 

7 = 7 1 + V m V m7 (2.4) 

where V m denotes derivation with respect to the slope dependence of 7 and 1 is 
the unity matrix. 

Strictly speaking this derivation applies only to equilibrium relaxation. For a 
far from equilibrium process like MBE growth it does not make sense to define a 
surface free energy per area 7(V/i), although the "true" contribution of the surface 
current due to relaxation by curvature most likely has a functional form of that 
type. The identification of these contributions and the correct derivation of the 
corresponding continuum form remains an open question. 
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2.3. Deposition, desorption and noise 

For most applications it will be correct to assume a deposition intensity F constant 
in time and space. Interesting cases with explicit space and time dependence are 

(i) pulsed molecular beams, used to decrease the second layer island nucleation 
rate, and (ii) geometric inhomogeneities if the beam does not reach every part 
of the surface under the same angle. The most important mechanism for implicit 
deposition heterogeneity is steering, where the incident atoms are deflected by 
attractive interaction with surface features such as islands and mounds, and effect 
which is of course most prominent under grazing incidence 01 I25| . 

At sufficiently high temperatures desorption may become an important ef- 
fect. Terrace adatoms then will be most susceptible to desorption at a rate 1/r, 
during their diffusion time before they reach the stronger binding steps. To a first 
approximation the desorption rate depends on the average distance to the next 
reachable step and therefore on the surface slope: the higher the slope, the more 
steps are present and the lower the desorption loss. Within an ansatz a la Burton- 
Cabrera-Frank .2! one can relate the desorption length x s — 2a\f~Dr, diffusion 
length £v, and slope V/i = a/£ to an expression for the desorption loss in 

volume per time and surface area 

£ W (x,f) = F [(x s /e)t & nh(£/x s ) - 1] « F y^^y?- (2-5) 

There are three sources for noise in MBE growth [21]: (i) Deposition is a 
stochastic process of single events, so F has to be complemented by a "shot noise" 
term £(x, i) with statistics 

(£(x, t)) = and (£(x, i)£(x', t')) = a 3 F <5(x-x') 5{t-t'). (2.6) 

(ii) Adatom diffusion also consists of single events, and to a first approximation 
surface currents are blurred by a "diffusion noise" term ^(x, t) with statistics 

(r/(x, t)) = and (ijfx, t) ■ r;(x', <')) = V 2 <5(x-x') 5(t-t'). (2.7) 

(iii) Nucleations also occur at random but the precise form of this contribution 
remains an open question. There are nontrivial correlations in time, as new islands 
tend to nucleate close to centres of previously nucleated islands. The space depen- 
dence of the nucleation noise correlator (f(x, i)z/(x', i')) is expected to have the 
form (V 2 ) 2 5(x-x') [33]. 

2.4. Structure and symmetry of continuum equations for MBE 

Heuristic surface currents as in the example in Fig. [21 lead to a symmetry in the 
resulting dynamical equation which can be used for evaluation, but is too strong 
a restriction compared to the "real" dynamics and may lead to erroneous results. 

If Jest uniquely depends on the slope m and not on derivative of its functions 
then crystal symmetry imposes JEST(m) = — Jest(— m). The resulting term of 
surface dynamics is invariant under the transformation h — > —h and so will be 
the surface configurations obtained from that. Breaking of this symmetry, e.g., by 
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terms as V/ ((V/i) 2 ) in J or an evaporation term e(((V/i) 2 ) do in fact change the 
characteristics of the dynamics substantially 1211 ITT] 

One can take advantage of this simple slope dependence if there is a scalar 
function V(m) such that Jest = — V m y(m). Stable zeroes of Jest then are 
minima of V(rxi). If additionally the curvature relaxation is of a simple "scalar" 
form a • 7 = Kl there is a Lyapunov functional of the dynamics 

~K 



2 -(Vm) 2 + ^(m) 



>2 



x (2. 



and the surface dynamics Eq. 1)2. 1J) can be written for the field of slopes m = V/i 
as 

d t m = d t {Vh) = VV ■— . (2.9) 
am 

If m is interpreted as a two-dimensional magnetic order parameter this from com- 
pares to conserved (model B) relaxational dynamics, with the only restriction that 
it remains a conservative field at all times, VAm=0. Facets in the surface are 
regions of constant slope, and appear as domains of constant "magnetisation" , 
edges between them are domain walls. As compared to its magnetic analogon, 
here the parallel component of m remains constant across a domain wall, which 
in the case of discrete stable zeroes of the current function (right panel of Fig. EJ) 
restricts the possible domain wall orientations and imposing constraints on the 
relaxational dynamics [1111191 12U] . 



3. Examples of applications 

Continuum equations are now applied in two examples of very basic choices in 
Eq. ([2.1[l . In some sense they highlight two opposite regimes, (i) on very long 
length and time scales coarsening of mounds in homoepitaxy, and (ii) on short 
scales resolving the lattice structure in the transition from layer-by-layer to rough 
growth. 

3.1. Coarsening 

In some sense the simplest example of a surface dynamics equation is the analogon 
to the classical XY-model in the sense of Eq. I[2.9II 

d t h = -{V 2 fh - V • [(1 - (V/i) 2 ) V/i] , (3.1) 

i.e., with the current function depicted in the left panel of Fig. [21 An initially 
flat surface h = const is linearly unstable to fluctuations of wave lengths A > 2tt 
which therefore initially grow exponentially. Once slopes of 0(1) are reached the 
surface organises itself into a pattern of regular mounds and troughs which keep 
the symmetry h <-> —h of Eq. (|3.1I) . The in-plane isotropy of (|3 - 1 f) is locally broken 
by the need to arrange the mounds and troughs into a regular pattern and facets 
with sharp wedges form |14| . 

Subsequently the pattern of mounds and troughs coarsens keeping a single 
typical length scale and a statistically self similar pattern at every time. The 
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Figure 4. Coarsening of surface mounds, between subsequent 
panels time advances by a factor of 10. The local surface slope is 
represented by colours, facets appear as regions of uniform colour. 
Spontaneous facet formation, although no anisotropy is imposed, 
J EST (m) = (1 - |m| 2 ) m (c.f. left panel in Fig.EJ. 



average mound size is found to increase with time as L ~ i 1 / 3 [111 I20| . In a weak 
sense this can be proved as an upper bound for an appropriate time average of 
the length scales L observed up to a give time t [H] . Here only a quick handwaving 
argument is given following 

The spatial average of the interface width W = (h 2 ) (where (h) = 0) increases 
due to Eq. I|3.1(l as 

I d t W 2 = \ d t (h 2 ) = -((Ah) 2 ) + {(Vhf) - (IV^I 4 ). (3.2) 

By the sketch of Fig. 13.11 the contributions of these terms are estimated to be of 
order 1/L: The facets are flat and the curvature is restricted to narrow strips of 
with O(l) at the domain boundaries, which have a relative weight of 1/L. Also, 
here the slopes are smaller than the stable value |V/i| = 1 attained at the facets. 
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Figure 5. Sketch of facets and wedges, i.e., domains of constant 
slope and their separating boundaries, as obtained in the simula- 
tions shown in Fig. 12.41 The only macroscopic length scale is the 
domain size L. 



So all terms give a contribution of 0(1/ L) and because of the stable slope W ~ L 
the increase in width leads to the estimate 

d t L 2 ~ d t W 2 ~ 1/L, integrated to L~t 1/3 . (3.3) 

The same rough heuristic argument holds for anisotropic surface currents with 
discrete minima, i.e., a discrete set of preferential surface orientations, but the 
four-fold symmetry (right panel of Fig. EJ) plays a special role because two kinds 
of edges occur, such that in general the assumption of a single length scale L in 
the system is wrong 21) . 

3.2. Roughening 

Layer by layer growth is technologically interesting because one can monitor the 
subsequent deposition of layers by oscillations in the surface width. It is nearly 
flat around integer fillings and the roughness takes maxima at half integer fillings. 
During growth the interface roughens and the width oscillations decay [o|. The 
Ehrlich-Schwoebel barrier and mounding instability is one prominent mechanism 
destroying layer-by-layer growth, and so can be any inhomogeneity in the deposi- 
tion intensity F{ h y(x, t). Although weaker, even the fluctuations of shot noise have 
this effect, and it can be captured by a simple renormalisation group calculation. 

The simplest way to implement layer lattice effects into an continuum surface 
equation like Eq. (|2.1|l is pursued by the conserved sine-Gordon-equation 

d t hhc,t) = -KA 2 h( Xl t) - yAsin 27r/ ^ X '^ + F + £(x,t). (3.4) 
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Among other the surface current is driven to regions of incomplete fillings such 
that h(x,t)/a± is favoured to take integer values ^7j. The notation distinguishes 
between the lattice constant normal (a±) and parallel (ay) to the surface terraces 
which will be used in the following analysis. 

In a momentum shell renormalisation one can average over the fluctuations 
caused by the short wavelength components of the noise, £ = £ + <5£, where £ 
only contains modes of wavenumber smaller than some cutoff A and <5£ the others. 
One now tries to approximate the dynamic equation for the averaged height field 
h = (h)g£ by an equation like 13.4(1 with suitably adjusted parameters. After 
rescaling, an effective renormalised equation is obtained. Three important physical 
interpretations can be obtained from the renormalised parameters of Eq. (|3.4H . 

(i) by the interplay of the diving force F and the lattice potential the up-down 
symmetry (h «-» —h) is broken and a term — AA(V/i) 2 is created. 

(ii) The strength of the lattice potential V is damped out decays exponentially 
with coverage 8 — Ft like oc exp ^— a^/i^y/O^ 

(iii) In a finite size sample smaller than the so called layer coherence length £ = 
^Y)/\\ renormalisation has to stop before the lattice potential has decayed. Shot 
noise fluctuations cannot accumulate to the point where width oscillations are 
completely blurred. 

4. Discussion 

In presenting the model derivations and examples of their application it becomes 
clear that continuum models are a good tool for questions where the inherent sim- 
plifications are justified. The insight gained from partial differential equations is 
different from results of more microscopic model. Last not least continuum equa- 
tions can handle larger systems and longer times than Monte Carlo simulations or 
even Molecular Dynamics. 

In many cases the phenomena of particular technological interest cannot be 
captured by continuum equations of the type described in this lecture. Magic 
island shapes, surface reconstruction, step bunching, formation of quantum dots 
and wires are just a few examples. 

Even if one isn't interested in these inherently discrete phenomena them- 
selves they nevertheless carry over to the behaviour on larger scales and a naive 
straightforward derivation of large scale continuum models may be wrong. The im- 
plementation of noise for some questions, such as roughening, is one such example. 
Also, the effective "large scale" Ehrlich-Schwoebel barrier of a step is influenced by 
kinks along the step which may serve as channels for downward incorporation of a 
terrace adatom. Obtaining the correct effective barrier, dependent on orientation 
and curvature of the step certainly is a nontrivial task. 

An important part of the surface current J flows along steps, as has become 
clear in the above derivation. In ref. |15| the step current due to kink ES barriers 
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is derived as a function of the surface orientation or slope. An open question is the 
role of step curvature and therefore of surface curvature on the step current. 

Also on the surface itself, the regularising term as it is derived in Eq. 1)2. 
should be derived from true nonequilibrium arguments. To do so, one would have 
to incorporate island nucleations, the creation and annihilation of steps. One pos- 
sible approach is to introduce fields for the densities of adatoms, advacancies, 
steps, kinks, and for their orientations. Locally they could be treated as mean 
field densities with certain effective "reaction" rates, but in their spatial interac- 
tion, by diffusion and advective transport, one could obtain a conceptually more 
fundamental continuum model for surface growth. 

5. Conclusion 

In this lecture an overview of continuum models for MBE surface growth has been 
given which includes 

• the conceptual foundation of continuum models as a limit of more detailed 
small scale approaches, 

• heuristic derivations due to lack of rigorous or even exact ways to obtain 
them, 

• the structure and the symmetry of the equations obtained in that way, 

• two examples emphasising two opposite time limits in the growth process 
with the application of two contrary methods, 

• criticism of weak points in continuum modelling, open questions and some 
speculations about answers to them. 

The main purpose was to highlight the physical background and give some intuition 
of possible approaches to it as a basis for mathematical treatments, both rigorous 
and numerical. It should be seen as a complement of the lectures pQ and ^01 
which give a similar introduction to lattice and Burton-Cabrera-Frank models 
respectively. Although continuum models seem to fail in some important physical 
aspects they nevertheless are a good tool and subject for mathematical approaches. 
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